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Theory and simulations of toroidal and rod-like structures in single-molecule 
DNA condensation 
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DNA condensation by multivalent cations plays a crucial role in genome packaging in viruses and sperm 
heads, and has been extensively studied using single-molecule experimental methods. In those experiments, 
the values of the critical condensation forces have been used to estimate the amplitude of the attractive 
DNA-DNA interactions. Here, to describe these experiments, we developed an analytical model and a rigid 
body Langevin dynamics assay to investigate the behavior of a polymer with self-interactions, in the presence 
of a traction force applied at its extremities. We model self-interactions using a pairwise attractive potential, 
thereby treating the counterions implicitly. The analytical model allows to accurately predict the equilibrium 
structures of toroidal and rod-like condensed structures, and the dependence of the critical condensation force 
on the DNA length. We find that the critical condensation force depends strongly on the length of the DNA, 
and finite-size effects are important for molecules of length up to 10®/rm. Our Langevin dynamics simulations 
show that the force-extension behavior of the rod-like structures is very different from the toroidal ones, 
so that their presence in experiments should be easily detectable. In double-stranded DNA condensation 
experiments, the signature of the presence of rod-like structures was not unambiguously detected, suggesting 
that the polyamines used to condense DNA may protect it from bending sharply as needed in the rod-like 
structures. 
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I. INTRODUCTION 

Single-molecule DNA micromanipulations pro¬ 
vide a very powerful tool for the study of DNA 
mechanicsii^, DNA-DNA interactions^^— , and DNA- 
protein interactions^. In the case of DNA-DNA 
interactions, it was possible to study the reentrant 
behavior of cation-condensed DNA^i^i^, and measure 
accurately the free energy of DNA condensation^. 
These studies are of great importance to determine the 
behavior of condensed DNA, which is crucial for DNA 
packaging and genomic ejection from viruses^^, and also 
to provide an accurate estimate of the forces involved. 

In single-molecule DNA condensation experiments, a 
DNA molecule is tethered at one end to a surface, 
and at the other end to a micron-sized bead, which is 
trapped by magnetic or optical tweezers. The experi¬ 
mental extension-force curves show that there is a critical 
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traction force, above which the system is fully extended 
and behaves as a worm-like chain semiflexible polymer, 
and below which the system progressively folds to a com¬ 
pletely condensed state^^— di. This critical condensation 
force allows to estimate the attractive DNA-DNA inter¬ 
action per unit length. So far, this was done by neglecting 
finite size effects for long A DNA molecules. 

DNA condensation in single-molecule experiments was 
already studied theoreticallyi^d^, but the geometry of 
the condensed structure was always assumed a priori to 
be toroidal. It is possible however that toroids are not the 
only possible geometry for the condensate, as was already 
noticed in theoretical studies of condensation of semiflex¬ 
ible polymers with self-interaction o^^d^ . It is therefore 
interesting to explore the possibility of the appearance of 
other geometrical configurations in single-molecule DNA 
condensation experiments. 

In this paper we propose to study DNA condensa¬ 
tion in single-molecule experiments by means of Langevin 
dynamics (LD) simulations, together with an analytical 
model. The scope of this work is two-fold: on the one 
hand, to better account for finite-size effects in the es¬ 
timate of attractive DNA-DNA interactions using mea¬ 
sured critical condensation forces; on the other hand, to 
evaluate the possibility of the appearance of the rod-like 
geometry of the condensate, alongside the well-known 
toroidal one. 

Our LD simulation method is based on the work 
of Carrivain et ali^, which has been successfully ap¬ 
plied to modeling single-molecule DNA supercoiling 
experiments^. The main advantage of using this tech¬ 
nique is its speed and enhanced sampling efficiency, based 
on (a) efficient algorithms to simulate rigid bodiesi^ and 
(b) the global thermostat scheme introduced by Bussi 
and ParrinelloiS. Using this method, we are able to ob¬ 
serve nucleation of DNA condensation within our Simula- 
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tion windows. The well-known toroidal structure of DNA 
condensates^ is not the only geometry that we observe: 
rod-like structures also appear. 

The analytical model we propose is based on the the¬ 
oretical ansatz proposed by Hoang et al3i. This was 
developed to study toroidal and rod-like structures of 
free DNA in solution, during crowding- or cation-induced 
condensation. We extended the model proposed by these 
authors to include the presence of a traction force ap¬ 
plied at the extremities of the DNA. Using this model, 
we are able to evaluate the different contributions to the 
free energy of the system: the bending energy, the sur¬ 
face tension, the bulk attractive force, and the worm-like 
chain entropic free energy of the non-condensed region 
of the DNA. We are also able to calculate the critical 
condensation force at any DNA length. 

In both our simulations and our analytical model we 
consistently model the counterions that are present in 
solution implicitly, using a pairwise attractive potential. 
The validity of this approximation is supported by a fa¬ 
vorable comparison of our results with experimental data 
on DNA condensation. 


II. METHODS 
A. Geometry 

The geometry of our system is described in figure [T] 
We describe a polymer of length L, tethered at one end 
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FIG. 1. Geometry of the single-molecule DNA manipulation 
experiment.(A) Bead and DNA molecule attached to it; (B) 
rod-like condensed structure: R is the rod diameter, 7 is the 
aspect ratio of the rod; (G) toroid-like condensed structure; 
R is the toroid radius, a is the ratio between the radius and 
the thickness of the toroid. 

to a surface and at another end to a bead. In magnetic 
tweezer experiments, the force F acting on the bead is im¬ 
posed, and the end-to-end distance z is measured. In op¬ 


tical tweezers, the bead is trapped around a fixed height 
zq, and the force is measured. These two experimental 
setups correspond to different statistical ensembles: at 
fixed force (where we use the Gibbs free energy), and 
at fixed end-to-end extension (where the Helmholtz free 
energy is used). 


B. Analytical model 

We describe the polymer as a made of two different 
phases: a condensed phase characterized by a length 
and an extended phase of length L — Lc- We therefore 
have the following two expressions for the free energy, at 
fixed force and at fixed end-to-end extension respectively: 

G = Ec + {L - Lc)gwLc{F)j (1) 

A = Ec + {L — Lc)awLC ^ ^ , (2) 

where E^ is the energy of the condensed phase, gwLC 
and awLC are the free energies per unit length of the 
worm-like chain phase, in the fixed force and fixed exten¬ 
sion ensembles, respectively. In our model, we neglect 
the entropic contribution due to fluctuations inside the 
condensed phase: we consider that the condensed phase 
has a definite conformational structure. We now discuss 
separately the contributions from the condensed phase 
and from the extended phase. 

For the condensed phase, we base our model on the 
recent work by Hoang et ali^. Here, we report the main 
results of this study. More details on the derivation of 
the following formulas, as well as a discussion on their 
validity, is reported in the original study^i. 

The authors model a polymer chain with self¬ 
interactions as a chain of spheres of radius b connected 
by rigid bonds, and propose an energy ansatz for the 
toroidal and rod-like geometry. In our model, we take 
their expressions, and assign the energy of the condensed 
phase Ec to one of the two formulas: 

Etorozd{v,ct:Lc) = 2lpTr'^^^ + 

- </>(d)d-l5l/3^4/3^-l/3^-2/3^2/3^ 

+ 3(j3{d)b-^ Lc (3) 

Erodiv, 7, Lc) = Y^pT] [4r] (7 -k 4/3)]“^''^ 

- 2'n:(j){d)d~^b^/^ (7 -k 2 ) [Arj (7 -k 4/3)]“^^^ L^^^F 

-k S(j){d)b~^Lc- (4) 

Both the toroid and the rod geometries are character¬ 
ized by 77 , the packing fraction in the condensate. The 
toroid is described also by a, which is the ratio between 
the toroid radius and the thickness of the toroid; the 
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rod is described by 7 , the rod aspect ratio between the 
width and the length of the rod (see figure [T|). The inter¬ 
action of the spheres is given by the potential function 
4>{d). We chose the Lennard-Jones potential energy func¬ 
tion, which is numerically very close to the more realistic 
Morse potential: 


Hd) = e 


© 
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( 5 ) 


invert the pwLc{F) function and we obtain the FwLcip) 
function. Once equation IH] is specified for F = Fwlc{p), 
we finally obtain 

awLC = 9 wlc [FwLcip)] - FwLcip)p- ( 10 ) 

The equilibrium condition is obtained by minimizing 
the free energy (either equation [T] or [^l with respect to 
all the variables. 


Here, a is the Lennard-Jones radius of the monomers. 
The variable d, which represents the inter-monomer lat¬ 
eral distance in the condensed phase, is related to rj by 
the following expression: 


d — ^\/Phex / Pj 


( 6 ) 


We may derive interesting predictions from minimiza¬ 
tion of the Helmholtz free energy (equation [2]). The min¬ 
imum condition for is given by 


dE, 

dLc 


— gwLc 

Lc=L* 


Fwlc 



= 0 , 


( 11 ) 


where rjhex = 7 r/ 6 cot( 7 r/ 6 ). Finally, Ip is the polymer 
bending persistence length. Notice that the precise form 
of the interaction potential was found to be of minor 
importance in determining the position and length de¬ 
pendence of the toroid/rod transition line in the work 
of Hoang et al^i. We also checked that the Morse po¬ 
tential gave little or no difference compared to the the 
Lennard-Jones potential in our analysis. 

Both the toroid and the rod energy functions are a sum 
of three terms: the bending energy to fold the polymer 
into its condensed state (which is ~ Lc ); the surface 
tension (which is ~ and the bulk attractive en¬ 

ergy, which is proportional to The surface tension 
and the bulk energy terms are proportional to the am¬ 
plitude of the attractive potential, e. These terms are 
derived by supposing that the packing of the polymer in 
the condensate is hexagonal, so that one can calculate 
the average number of neighbors in the bulk, and the av¬ 
erage number of molecules at the surface. As we shall 
see, the dependence of these terms on has important 
consequences on the equilibrium properties of the system. 

The free energy of the extended phase may be ex¬ 
pressed using the formula derived by Marko and Siggiai^^ 
for the Gibbs free energy per unit length of a worm-like 
chain semi-flexible polymer of length L — Lc: 

) = mta { (^ - f) (cotl.2a - T) | . (7) 

This expression is used to calculate the relative extension 
p = z/{L — Lc) at thermal equilibrium 


(cnx dgwLC 
Pwlc{F) = -. 


( 8 ) 


To express the free energy per unit length in the fixed 
end-to-end extension ensemble, we must perform the Leg¬ 
endre transform of gwLC- 


where the starred symbols refer to the equilibrium con¬ 
dition. The above equation is derived in the Appendix. 
In the thermodynamic limit, Lc —> oo, and equation 1111 
may be approximated by 

- 3 £ = gwLC [Fwlc ( p *)] • ( 12 ) 

This equation depends only on p* = z/{L — L'^), and 
not on z and L separately. Therefore, in the limit of 
large L we expect z/{L — L*) to be constant. This may 
be translated into the following useful formula, which we 
will use later: 


L 


★ 

C 



(13) 


At fixed end-to-end extension, the force acting on the 
bead at thermal equilibrium is given by: 


F = Fwlc [p*] • (14) 


This equation is derived in the Appendix. 

To obtain the equilibrium structures, we use numerical 
routines for minimization of multi-variable functions, as 
implemented in the GNU Scientific Library^^. 

To calculate the critical condensation force, we switch 
back to the fixed force ensemble. The transition from 
the fully extended state to the fully condensed state is 
obtained by equating the chemical potentials of these two 
phases, at the critical force. We then obtain the following 
expression: 


= 9wlc{Fc), (15) 

where in E* we set Lc = L. This equation may be solved 
numerically to obtain Fc- Since the bending energy and 
surface tension terms vary as ~ lJ^^ and ~ Lc^^ respec¬ 
tively, in the thermodynamic limit of L —>■ oo, we obtain 


dwLC = 9 wlc{F) - Fp. (9) 


- 3e = 9 wlc{F^), (16) 


To obtain the force as a function of the extension, we 


where F^ is the critical force for an infinitely long chain. 
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C. Molecular dynamics simulations 

Our simulation assay is based on the recently published 
work by Carrivain et ali^. We mention here only the most 
relevant features of this simulation method. First, it is 
a simulation of rigid bodies, where the angular velocity 
is taken into account, and the simulated bodies have an 
impenetrable volume. In our case, the DNA is composed 
by N cylinders of length Ig, with a hard-core radius equal 
to the DNA crystallographic radius, and is attached to a 
magnetical or optical bead, which is treated as a sphere 
(see figure [T]4) . Second, the DNA chain is modelled by 
an articulated system, where the joints act as holonomic 
contraints of a mechanical system. Forces and torques 
may be applied at each joint, so to model bending and 
twisting energies. In our case, we apply only a bending 
term, and neglect the twisting contribution, which is ir¬ 
relevant because the DNA molecules are typically nicked 
in the experiments. Finally, the coupling to the ther¬ 
mal bath is treated using the global thermostat scheme, 
introduced by Bussi and ParrinelloiS. This allows for 
fast equilibration of the system and realistic treatment 
of thermal fluctuations. 

We model the bending energy associated to two suc¬ 
cessive cylinders as: 

Eb{0) = (17) 

where 0 is the angle formed by the tangent vectors of the 
two cylinders, and gb is a constant related to the bending 
persistence length (for full details, see referencei^) . 

Here, we include the effect of intra-DNA interactions, 
as modelled by the Lennard-Jones potential equation 0 
In our system, the potential acts between the centers 
of mass of the cylinders that compose the articulated 
system. The Lennard-Jones radius a was chosen so that 
the minimum of the potential corresponds to 28A. The 
full list of parameters used in our simulations is given in 
table m 

To simulate optical tweezer experiments, we add a po¬ 
tential energy term to the bead, expressed as an isotropic 
harmonic trap: 

Etrap — '^^trap{\‘^ ^o|) ■ (1^) 

Here, ktrap is the stiffness of the optical trap. Typical 
values of this parameter are of the order of O.lpN/nm^. 
The vectors r and Tq = (xo,yo, Zq) are the actual position 
of the bead, and the position of the center of the force 
(which would correspond to the laser focus), respectively. 

We run an equilibration round of 10® steps, using a 
local Langevin thermostat and no inter-monomer poten¬ 
tial. After this, we switch to the global thermostat and 
turn on the inter-monomer interactions, and perform 10^ 
LD steps for production run. 


D. Data analysis 

We briefly describe here a few methods we used to 
analyze the data from our LD simulations. 

Condensate length. To calculate the length of the 
condensed region from an LD trajectory, we proceed as 
follows. First, we calculate the distance matrix between 
all the DNA segments at a given time step. That is, 
calculate — rj\. Then we select all segments in 

which \i—j\ > 1, and Mij < Ig. The number of segments 
that satisfy this condition, multiplied by h, gives Lc- 
Free energy. The free energy of a given state may 
be estimated by calculating the total energy of the con¬ 
densate (Lennard-Jones potential plus bending energy), 
and adding the worm-like chain contribution of the un¬ 
condensed region. This can be calculated by averaging 
Lc over the trajectory, and using the expression for the 
worm-like chain free energy per unit length proposed by 
Marko and Siggia (see equations [TUI and 171). As we have 
done in our analytical model, we neglect the entropic con¬ 
tribution of fluctuations of the segments in the condensed 
phase. 

Critical force. It is challenging to estimate Fc with LD 
simulations directly. At the critical point, the time scale 
for the formation of an initial condensation loop diverges. 
Therefore, it is difficult to estimate Fc by starting from a 
high force value and decreasing it, because the time scale 
for nucleation is inaccessible in LD simulations. On the 
other hand, a similar problem is encountered when start¬ 
ing from a condensed structure and increasing progres¬ 
sively the applied force. Here, the kinetic barrier between 
the folded and unfolded states is significant, therefore 
presenting the same conceptual problem. However, one 
can estimate Fc by fitting the (Lc) {{z)) curves to a line, 
and estimating the intercept of the line with the Lc = 0 
axis (see equation m and figure El). The value of z at 
which Lc = 0 gives an estimate of p*, from which we can 
estimate the critical force using the Fy/Lc{p) function. 
Distinguishing between toroids and rods. To au¬ 
tomatically detect whether a condensed state is toroidal 
or rod-like, we developed an heuristic method outlined 
below. Once the cylinders that are part of the condensed 
phase are identified, for each pair of cylinders i and j in 
contact, we evaluate cosyij = F ■ tj. If cosjij > 0, we 
say that those segments are in parallel contact, otherwise 
they are in antiparallel contact. We can then define an 
order parameter as s = {Np — Na)/Nc, where Np and 
Na are the number of parallel and antiparallel contacts, 
respectively, and Nc is the total number of contacts. For 
a toroid, most segments will be in parallel contact, so we 
expect s w I. On the other hand, for a rod-like configu¬ 
ration we expect s r; 0, since Np ss Na- 

III. RESULTS 

We performed Langevin dynamics simulations of a 3 
kb DNA molecule {N = 150, L = Ipm), as described in 
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Methods, using parameters shown in tablelH in an optical 
tweezer-like setup. We compared the simulation results 
with the predictions of our analytical model (equations 
[Hand 121) using the same parameters, for the toroidal and 
rod-like geometries. We present here the results of these 
two approaches. 

A. Fixed zi^jL simulations 

First, we performed 20 independent simulations of the 
DNA chain, at fixed zq/L = 0.6 and at e = O.TfcsT/nm. 
We observed that after 10^ time steps, 9 chains adopted 
the rod-like condensed state, 10 adopted the toroidal 
state, and one remained uncondensed. Snapshots of the 
rod-like and toroidal states are shown in figure [TJl and 

dp. 

Figure [2] shows two example traces of {L — Lc{t))/L for 
a toroid and a rod. The condensation of the chain occurs 
as discrete jumps, both for the toroid and rod-like ge¬ 
ometries. The steps correspond to adhesion of successive 
portions of the DNA chain onto an initially formed con¬ 
densation loop for toroids and a condensation stretch for 
rods. In the time trace of the toroid condensation, a more 
progressive, linear condensation step is also observed. 



FIG. 2. Variation of the uncondensed length {L — Lc)/L {Lc is 
calculated as explained in Methods, as a function of time for 
two specific simulations, one which nucleates into a rod-like 
condensate (solid red) and one which nucleates into a toroid 
(dashed blue). 


B. Force-extension curves 

Next, we used the final configurations of the simula¬ 
tions described in the previous section to perform stretch¬ 
ing simulations. Here, every 10® time steps the value 
of zq/L was incremented by 5%. The resulting force- 
extension curves are depicted in figure [31 The toroid- and 




FIG. 3. Average force versus average extension of single 
stretched rods (A) and toroids (B). The different symbols 
correspond to simulations starting from independent config¬ 
urations (see Main Text). The worm-like chain behavior for 
a completely uncondensed DNA is shown as a black continu¬ 
ous line. In A, black dashed lines correspond to the fit to a 
worm-like chain behavior for, from left to right, L/la = 105.5, 
L/ls = 108.5, 118.5, 122.5, 132.5. The simulation protocol is 
described in Methods, and the parameters used are shown in 
table m 

rod-like structures have a markedly different response un¬ 
der stretching force. The toroidal structures unfold pro¬ 
gressively, exhibiting a saw-tooth-like force-extension be¬ 
havior, eventually becoming fully extended (see final blue 
points in figure|3]). The peaks in the force value can go up 
to ^30pN. In contrast, rod-like structures never unfold, 
even at the highest forces probed here (^^SOpN). In all 
our simulations, we could not unfold rod-like structures, 
whereas toroidal structures all eventually unfold. 

The difference in the force-extension behavior of 
toroids and rods under tension can be understood in 
terms of the direction of the applied force. In a rod-like 
condensate, the rod axis is in general approximately par- 
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Parameter Description Value 


Ip 

Bending persistence length 

45 nm 

T 

Absolute temperature 

300 K 

N 

Number of DNA segments 

150 

Us 

Number of base pairs per DNA segment 20 

n 

Thermostat coupling frequency 

3 10i°s-i 

ktrap 

Optical trap stiffness 

0.2 pN nm 

e 

Lennard-Jones potential well depth 

0.2-0.7 ksT/nm 

a 

Lennard-Jones radius 

25.08 A 


TABLE I. Summary of the parameters used in Langevin dynamics simulations and in our analytical model (where applicable). 


allel to the direction of the applied force (see figure [TJ3). 
As the force increases, the direction of the force becomes 
more and more orthogonal to the direction of the DNA- 
DNA contacts. This makes it difficult for the force to 
unfold the condensate. On the other hand, the toroidal 
structures have a larger degree of lateral mobility, that 
results in fluctuations aligning the lateral contacts with 
the direction of the applied force. 


C. Free energy of toroidal and rod-like structures 

To further gain insight on the system, and to test the 
results of our analytical model (equations [T] and [2]) , we 
performed simulations at different values of zq/L. For 
each value of zq/L, we obtained results for 10 indepen¬ 
dent starting configurations. We then calculated the free 
energy of the final configurations as outlined in Meth¬ 
ods, averaging over the 10 simulations. Figure U] shows 
the comparison of the numerical minimization of our an¬ 
alytical model, for the toroid and rod geometries, and 
the free energies obtained from Langevin dynamics sim¬ 
ulations as outlined above. 

For each value of zo/L that we studied, the toroidal 
states are the ones having the lowest free energy. The free 
energy minima for toroidal and rod-like configurations, 
however, are remarkably close (within ~ Q.lksT/mn dif¬ 
ference). This fact is also manifest in the observation of 
simulations in which a toroidal and a rod-like condensed 
phase coexist within the same DNA chain (not shown). 

Our simulations show that at {z) jL > 0.7 for toroids 
and {z) jL > 0.6 for rods, the condensed state ceases to 
exist. At values of {z) jL higher than this, the free energy 
is equal to the worm-like chain free energy. 

We notice that there is a significant spread between 
the values of the free energies calculated in our simula¬ 
tions. However, the lowest values of the free energy for 
a given conformation are remarkably close to the pre¬ 
dictions of our analytical model. This means that in the 
other cases, the conformation adopted by the chain is not 
an equilibrium one, but a local minimum. This is a sign 
that the final conformation is strongly dependent on the 
nucleation of a first condensation loop. 



FIG. 4. Minimum free energy of rod-like and toroidal geome¬ 
tries (blue and red dashed lines, respectively), as calculated 
from numerical minimization of equation [2l and average free 
energies of final states in Langevin dynamics simulations, cal¬ 
culated as explained in Methods, for rods (red circles) and 
toroids (blue triangles). 


D. Critical condensation force 

Figure [5] shows the values of the critical condensa¬ 
tion force Fc as estimated from the analytical model 
(see equation [15]), as a function of the DNA length, at 
fixed e. The two values chosen for e correspond to the 
case of spermidine (0.20 kETjum) and spermine (0.33 
ksT /nm) given in reference^. As we explained in the 
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Methods section, for a very long chain the critical con¬ 
densation force reaches a finite value , and becomes 
independent of the geometrical structure of the conden¬ 
sate. Here, we show that this thermodynamic limit is 
reached only for very long molecules. As an example, 
for e = 0.33fcsT/nm, the value of Fc reaches F^ only 
for L = 10®nm, the size of a genomic DNA molecule. 
The experimental values of the maximum condensation 
force, as extracted from the work of Todd et al^ and Mu- 
rayama et al^ are also depicted. There is good agreement 
between the theoretical and experimental values. 



, ^ 

10 ° 10 ^ 10 ^ 10 ° 10 ^ 10 ^ 
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FIG. 5. Dependence of the critical condensation force Fc 
on the length of the DNA, for e = O.SSfcsT/nm (A) and 
e = 0.20fesr/nm (B). These two values correspond to the 
case of spermidine and spermine, respectively. The points are 
extracted from the experimental data of Murayama et al^, 
Besteman et al^, and Todd et al^. 


linear dependence of < Lc > on < z >, which is one of the 
key predictions of our analytical model (see Methods). 
A linear fit of this data provides a value of the intercept 
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FIG. 6. Gondensed length in simulations, as a function of 
average (z) jL. Condensed length was calculated as described 
in Methods. Error bars are standard deviations. The dashed 
lines are linear fits to the data. 

with the z = 0 axis which is close to I, meaning that in 
the absence of an external force, the whole molecule is in 
the condensed state. The value of p* estimated from our 
data is 0.96 and 0.87 for toroids and rods, respectively. 
Unfortunately, small errors in p* result in large errors on 
the estimate of F^ as for p > 0.85 the Fwlc{p) function 
has a very steep increase (see figure [3] for comparison). 
Therefore, we estimate F^ to be between Ril.5 and 20pN, 
at £ = 0.7A:Br/nm. 


The reason for the very slow convergence of Fc to F^ 
is to be found in the surface tension term. In fact, this 
term is proportional to e (see equations [3] and |H which is 
typically a very significant value, and depends on Lp , 
which is not far from linear dependence. Therefore, this 
term represents a significant energy penalty even for a 
very long DNA. 

As outlined in the Methods section, one can estimate 
the critical condensation force from LD simulations by 
looking at the dependence of the condensed length Lc 
on z. Figure [B] shows the variation of the average Lc, 
for toroids and rods, at fixed z. The data clearly shows a 


IV. DISCUSSION 

In this study, we developed an analytical theory and a 
simulation method to investigate the behavior of a DNA 
under traction and in the presence of self-attraction. 
This study was primarily aimed at the investigation of 
single-molecule DNA condensation in optical and mag¬ 
netic tweezers^^— Whereas our simulations allow 
to study the kinetic aspects of DNA condensation, our 
analytical theory is aimed at the study of the equilibrium 
properties of the system. The latter are experimentally 

















accessible, using a very slow unloading rate of a mag¬ 
netic bead (see e. g. Refi^). In this section we discuss 
our results in relation to the experimental data. 


A. General features of single-molecule DNA condensation 
experiments are reproduced by our simulations 

Firstly, we notice that several experimental features 
of DNA condensation are well captured by our model. 
The condensation of the chain proceeds in discrete steps, 
and each step was shown to correspond to folding of 
a portion of the molecule onto an initial condensation 
loop^i^i. When the condensate is toroidal, the force- 
extension curve has the characteristic saw-tooth shape, 
which was seen in experiments^ and in previous simu¬ 
lations of DNA condensationi^. The force peaks of our 
simulations are quantitatively in agreement with the ex¬ 
perimental ones, being in the 10-30 pN range (see figure 
[SI). Our values are somewhat higher than the experimen¬ 
tal ones, and this may due to the limited time window 
accessible by LD simulations. Moreover, we find that 
the conformation of the condensed structures apprearing 
in the LD simulations is strongly dependent on the size 
of the initial condensation loop, which is a well-known 
feature of DNA toroids observed in electron microscopy 
experiments2^“— . 

Despite its simplicity, our simple assumption of a pair¬ 
wise potential between the monomers is able to capture 
some of the principal features of the system. 


B. Evaluating the amplitude of DNA-DNA interactions 
using values of critical condensation force 

The analysis of the dependence of the critical conden¬ 
sation force on the length of the DNA gives an interesting 
prediction: the limit of infinitely long chains is reached 
very slowly, as shown in figure |S] In the work of Todd et 
al^i^, the authors accurately measured the condensation 
force at different ionic concentrations. The condensation 
force has a maximum value as a function of the concen¬ 
tration of multivalent ions. For this value, one can show 
that the entropic contribution to the free energy due to 
ion mixing vanishes, and one can estimate the amplitude 
of the attractive DNA-DNA interactions using the value 
of the maximum critical force. This estimate is based on 
the idea that for a very long DNA molecule, the total 
free energy of the condensate may be approximated by 
its bulk value. Our study provides a quantitative way of 
veryfing this hypothesis. Figure |3 shows the relationship 
between the estimated values of e and the condensation 
force. The dashed line corresponds to the hypothesis of 
Todd et al, i. e. assuming direct proportionality be¬ 
tween e and Fc, independently of the molecular length. 
Remarkably, the curve for L = 16.3^m (corresponding to 
the A DNA molecules used in the experiments of Todd 
et al) falls almost exactly on the dashed line. There¬ 


fore, even though e should have been estimated more 
rigorously, the values supplied by Todd et al for e are 
correct. Using the fact that those values are correct, we 



FIG. 7. Estimated value of the amplitude of the DNA-DNA 
attraction e, from a given measured value of the maximum 
critical condensation force Fc. Solid colored lines correspond 
to the theoretical prediction of our analytical model, as ex¬ 
tracted by numerically solving eauation llSI Dotted line corre¬ 
sponds to the case of L = oo. Dashed line corresponds to the 
estimate done by Todd et al^, assuming direct proportional¬ 
ity between e and Fc- Experimental points were assigned by 
assuming that the values of e estimated by Todd et al^(blue 
triangles) were correct (see main text for further details) at a 
given ionic condition, and assigning the points extracted from 
the studies of Murayama et al^ (green circle) and by Besteman 
et al2(red squares) to the corresponding ionic conditions. 

also set other values of experimental critical forces in the 
figure, as extracted from the works of Murayama et al^ 
and Besteman et al^, along with the theoretical estimates 
corresponding to the molecular lengths used in those ex¬ 
periments. The theoretical and experimental values are 
in reasonable agreement. 

C. Force-extension curves allow to discriminate between 
toroids and rods 

This study also showed that rod-like condensates have 
a different force-extension behavior compared to toroids. 
In the paper by Baumann et alJ^, it was reported that 
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some of the plasmid-length DNA molecules could not be 
stretched to the full original extension. Our results may 
provide a possible explanation for the two different force- 
extension behaviors seen experimentally. To our knowl¬ 
edge, the results of Baumann et al. are the only ones that 
show two different force-extension curve classes. In other 
experimental works^ on the same system, this was not 
reported. 

The geometrical configuration of the condensate is 
highly dependent on the model chosen for the DNA 
bending. If one allows for formation of kinks in the 
chain, with low energetic cost, the rod-like structure be¬ 
comes the most favorable geometry for the condensate. 
In contrast, adding a high energy cost for the hairpin 
turns strongly favors the toroidal geometry (data not 
shown). These results are in line with previous theo¬ 
retical studiesi^i^. This, together with the fact that the 
force-extension signature of rods was not clearly reported 
in the literature, suggests that the high bending angles 
required in the rod-like structures have a higher energetic 
cost than that expected by a simple harmonic model as 
we used here (see equation HZl). One possibility is that 
the polyamines used in DNA condensation experiments 
protect the secondary structure of DNA from kinking or 
bending sharply. It is in fact known that the polyamines 
stabilize base-pairing22i and do not disrupt the B-DNA 
structural fornt^. It is therefore unlikely that DNA kinks 
are favored by inhomogeneous patterns of adsorption of 
counterions on DNA. 

The structure of condensed DNA is a subject of great 
intereslji^, because of the importance, for instance, of 
DNA ejection from viral capsids^. Our results show 
that under the influence of an external traction force, 
the toroidal condensed structures are the only ones that 
can readily unfold. When a DNA molecule becomes rod¬ 
like, the unfolding is strongly impeded. This could be an 
important element to take into account when envisaging 
therapeutic strategies against DNA viruses. 
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Appendix: Derivation of equations 1111 and 1141 


The force acting on the bead at thermal equilibrium is 
obtained by taking the total derivative of the minimum 


Helmholtz free energy with respect to 


_ dA* ^ dEc 


dz 


+ 


Xi=X* 


dz 


{L — Lc)awLC 


L-L, 


(A.l) 


Lc=Lt 


where the starred symbols indicate the equilibrium values 
of the variables. In the above expression, the first term on 
the right-hand side vanishes, because the gradient of the 
condensed energy is zero at equilibrium. The second term 
is the only one that contributes to the force, and may be 
calculated by straightforward differentiation. The first 
step is to calculate the derivative of awLC with respect 
to p: 


d d 

-^awLcip) = ^ {gwLC [FwLcip)] + FwLc{p)p} = 


dPwLC dgwLC 
dp dFwLG 


dFwLC , p 
■ P —5 - 1 - Fwlckp) = 


dp 


= FwLcip), (A-2) 


where we have taken into account equation |8l Using this 
result we may now calculate easily the force acting on 
the molecule at equilibrium by inserting equation IA.2I 
into eauation lA.il 


= FwLciP*), 
(A.3) 

which is equation [TT] 

We may now explicitly calculate the derivative of the 
Helmholtz free energy equation [5] with respect to Lc- 




{L — Lc) awLC 


L-F 


dE, 

dLc 


Lr=L* 


d 


{L — Lc) awLC 


L-L, 


= 0 . 

(A.4) 


The second term on the left-hand side of this equation 
can be calculated using again equation lA.21 


d 

dLc 


{L — Lc) awLC 


L-Lc 


= -awLcip) + {L - Lc)i:^ir-awLcip) = 

C/J-JQ ^ P 

= -awLcip) + {L - Lc)— - i-^Fwlc (p) = 


= —gwLC 


{L - Lcf 


Fwlc 


^L-Lc 

Using this result, equation IA.4I becomes 


dFc 

dLc 


- gwLC [Fwlc (p*)] = 0, 


. (A.5) 


(A.6) 


F=L* 

which is equation [TT] 
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